clc; clear; close all;

Master-Thesis - Model and Calibrations:

Utility CASE 1

The utility functions are structured as follows:
U_ij = R + theta * s_i - p_ij - tau * (x - x_j)^2
where:
The four products are:
% The utility functions are defined as follows;
 
 
syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda
 
U_H1 = R + theta*s_H - p_H1 - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0 - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0 - tau*x^2;
 
disp(U_H1);
disp(U_H0);
disp(U_L1);
disp(U_L0);

Probabilities CASE 1

This section defines the choice probabilities for each product using a logit model with a scale parameter lambda.
The probabilities are derived from the utility functions of each product. The general form of the logit probability for product ij is:
Pr_ij = exp(lambda * U_ij) / Sum(exp(lambda * U))
Lambda is a scale parameter that controls the sensitivity of choice probabilities to utility differences. The denominator is the sum of exponentiated utilities for all choices, ensuring that the probabilities sum to 1. These probabilities represent the likelihood that a given consumer, characterized by their location (x) and willingness to pay (θ), will choose each product.
denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
 
disp(Pr_H1)
disp(Pr_H0)
disp(Pr_L1)
disp(Pr_L0)
 
% Pr_H1 = exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)));
% Pr_H0 = exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))/(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)));
% Pr_L1 = exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)));
% Pr_L0 = exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))/(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)));

Demand functions CASE 1

Demand defined as:
D_ij = int_0^1* int_theta_min^theta_max*Pr_ij(x, theta) f(x) g(theat) dtheta dx
where:
since it is integrated over all consumers the total demand or market share is captured by the intergation.
 
syms x theta p_H1 p_H0 p_L1 p_L0 theta_min theta_max
 
f_x = 1;
g_theta = 1 / (theta_max - theta_min);
 
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
 
disp('D_H1 = '); disp(D_H1)
D_H1 =
disp('D_H0 = '); disp(D_H0)
D_H0 =
disp('D_L1 = '); disp(D_L1)
D_L1 =
disp('D_L0 = '); disp(D_L0)
D_L0 =
 
% D_H1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% D_H0 = int(int(-exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% D_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% D_L0 = int(int(-exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);

Profit functions CASE 1

Profit for each product is defined as:
Pi_ij = (p_ij - c_ij) * D_ij
 
syms c_H1 c_H0 c_L1 c_L0
 
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
 
disp('Pi_H1 = '); disp(Pi_H1)
Pi_H1 =
disp('Pi_H0 = '); disp(Pi_H0)
Pi_H0 =
disp('Pi_L1 = '); disp(Pi_L1)
Pi_L1 =
disp('Pi_L0 = '); disp(Pi_L0)
Pi_L0 =
 
% Pi_H1 = -(c_H1 - p_H1)*int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% Pi_H0 = -int(int(-exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)*(c_H0 - p_H0);
% Pi_L1 = -(c_L1 - p_L1)*int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% Pi_L0 = -int(int(-exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)*(c_L0 - p_L0);

Profit maximization CASE 1

To find the profit-maximizing prices, we take the derivative of each profit function with respect to its price and set it equal to zero:
d(Pi_ij) / d(p_ij) = 0
dPi_H1_dpH1 = diff(Pi_H1, p_H1) == 0;
dPi_H0_dpH0 = diff(Pi_H0, p_H0) == 0;
dPi_L1_dpL1 = diff(Pi_L1, p_L1) == 0;
dPi_L0_dpL0 = diff(Pi_L0, p_L0) == 0;
 
disp('dPi_H1_dpH1:'); disp(dPi_H1_dpH1)
dPi_H1_dpH1:
disp('dPi_H0_dpH0:'); disp(dPi_H0_dpH0)
dPi_H0_dpH0:
disp('dPi_L1_dpL1:'); disp(dPi_L1_dpL1)
dPi_L1_dpL1:
disp('dPi_L0_dpL0:'); disp(dPi_L0_dpL0)
dPi_L0_dpL0:

Solve FOC for optimal prices CASE 1

solve the system of equations to optain optimal prices p_ij^*
% dPi_H1/dp_H1 = 0:
 
% int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/
% ((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))
% + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))
% + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)
% - int(int((lambda*exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))))
% - (lambda*exp(2*lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))^2)
% , theta, theta_min, theta_max), x, 0, 1)*(c_H1 - p_H1) == 0
 
% dPi_H0/dp_H0 = 0:
 
% int(int(-exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))/
% ((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))
% + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))
% + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)
% - (c_H0 - p_H0)*int(int((lambda*exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))))
% - (lambda*exp(2*lambda*(- tau*x^2 + R - p_H0 + s_H*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))^2)
% , theta, theta_min, theta_max), x, 0, 1) == 0
 
% dPi_L1/dp_L1 = 0:
 
% int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/
% ((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))
% + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))
% + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)
% - int(int((lambda*exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))))
% - (lambda*exp(2*lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))^2)
% , theta, theta_min, theta_max), x, 0, 1)*(c_L1 - p_L1) == 0
 
% dPi_L0/dp_L0 = 0:
 
% int(int(-exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))/
% ((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))
% + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))
% + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)
% - (c_L0 - p_L0)*int(int((lambda*exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))))
% - (lambda*exp(2*lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))^2)
% , theta, theta_min, theta_max), x, 0, 1) == 0
 
% to find the profit-maximizing prices p_H1*, p_H0*, p_L1*, and p_L0*. => define parameteres substitute into the equations and solve
 
syms p_H1 p_H0 p_L1 p_L0 theta x
 
% % Define parameter values
R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
theta_min = 0.2;
theta_max = 2.2;
c_H1 = 0.35;
c_H0 = 0.35;
c_L1 = 0.1;
c_L0 = 0.1;
 
denominator = (theta_max - theta_min) * ( ...
exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + ...
exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + ...
exp(lambda*(-tau*x^2 + R - p_H0 + s_H*theta)) + ...
exp(lambda*(-tau*x^2 + R - p_L0 + s_L*theta)) );
 
% FOCs:
foc_H1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- int(int((lambda * exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_H1 - tau*(x - 1)^2 + s_H*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) * (c_H1 - p_H1) == 0;
 
foc_H0 = int(int(-exp(lambda*(-tau*x^2 + R - p_H0 + s_H*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_H0 - p_H0) * int(int((lambda * exp(lambda*(-tau*x^2 + R - p_H0 + s_H*theta))) / denominator ...
- (lambda * exp(2 * lambda * (-tau*x^2 + R - p_H0 + s_H*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
 
foc_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- int(int((lambda * exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_L1 - tau*(x - 1)^2 + s_L*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) * (c_L1 - p_L1) == 0;
 
foc_L0 = int(int(-exp(lambda*(-tau*x^2 + R - p_L0 + s_L*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_L0 - p_L0) * int(int((lambda * exp(lambda*(-tau*x^2 + R - p_L0 + s_L*theta))) / denominator ...
- (lambda * exp(2 * lambda * (-tau*x^2 + R - p_L0 + s_L*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
 
solution = vpasolve([foc_H1, foc_H0, foc_L1, foc_L0], [p_H1, p_H0, p_L1, p_L0]);
 
disp('Optimal Prices:');
Optimal Prices:
disp(['p_H1* = ', char(solution.p_H1)]);
p_H1* = 0.97818251806516813907110255397178
disp(['p_H0* = ', char(solution.p_H0)]);
p_H0* = 0.97818251806516813907110255397178
disp(['p_L1* = ', char(solution.p_L1)]);
p_L1* = 0.67299421682353783913406045511087
disp(['p_L0* = ', char(solution.p_L0)]);
p_L0* = 0.67299421682353783913406045511087
 

Numerical results - CASE 1:

R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
theta_min = 0.2;
theta_max = 2.2;
c_H1 = 0.35;
c_H0 = 0.35;
c_L1 = 0.1;
c_L0 = 0.1;
 
p_H1 = 0.9781;
p_H0 = 0.9781;
p_L1 = 0.6729;
p_L0 = 0.6729;
 
 
denominator = @(x, theta) (theta_max - theta_min) .* ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) + ...
exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) );
 
D_H1 = integral2(@(theta, x) exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
 
fprintf('Demands:\n');
Demands:
fprintf('D_H1 = %.4f\n', D_H1);
D_H1 = 0.3166
fprintf('D_H0 = %.4f\n', D_H0);
D_H0 = 0.3166
fprintf('D_L1 = %.4f\n', D_L1);
D_L1 = 0.1834
fprintf('D_L0 = %.4f\n', D_L0);
D_L0 = 0.1834
 
fprintf('\nProfits:\n');
Profits:
fprintf('Pi_H1 = %.4f\n', Pi_H1);
Pi_H1 = 0.1989
fprintf('Pi_H0 = %.4f\n', Pi_H0);
Pi_H0 = 0.1989
fprintf('Pi_L1 = %.4f\n', Pi_L1);
Pi_L1 = 0.1051
fprintf('Pi_L0 = %.4f\n', Pi_L0);
Pi_L0 = 0.1051
 

Probability visualization for consumer choices - CASE 1

N = 10000;
R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
p_H1 = 0.97; p_H0 = 0.97; p_L1 = 0.67; p_L0 = 0.67;
 
x_samples = rand(N, 1);
theta_samples = 0.2 + (2.2 - 0.2) * rand(N,1);
U_H1 = R + s_H * theta_samples - p_H1 - tau * (1 - x_samples).^2;
U_H0 = R + s_H * theta_samples - p_H0 - tau * x_samples.^2;
U_L1 = R + s_L * theta_samples - p_L1 - tau * (1 - x_samples).^2;
U_L0 = R + s_L * theta_samples - p_L0 - tau * x_samples.^2;
 
denominator = exp(lambda * U_H1) + exp(lambda * U_H0) + exp(lambda * U_L1) + exp(lambda * U_L0);
 
P_H1 = exp(lambda * U_H1) ./ denominator;
P_H0 = exp(lambda * U_H0) ./ denominator;
P_L1 = exp(lambda * U_L1) ./ denominator;
P_L0 = exp(lambda * U_L0) ./ denominator;
 
figure;
subplot(2, 2, 2);
scatter(x_samples, theta_samples, 10, P_H1, 'filled');
title('Probability for H1');
xlabel('x');
ylabel('\theta');
colorbar;
colormap jet;
subplot(2, 2, 1);
scatter(x_samples, theta_samples, 10, P_H0, 'filled');
title('Probability for H0');
xlabel('x');
ylabel('\theta');
colorbar;
colormap jet;
subplot(2, 2, 4);
scatter(x_samples, theta_samples, 10, P_L1, 'filled');
title('Probability for L1');
xlabel('x');
ylabel('\theta');
colorbar;
colormap jet;
subplot(2, 2, 3);
scatter(x_samples, theta_samples, 10, P_L0, 'filled');
title('Probability for L0');
xlabel('x');
ylabel('\theta');
colorbar;
colormap jet;
 
% end of CASE 1:
clc; clear; close all;

Master-Thesis - Model adjusted for Tax on ICEV and Calibrations:

Utility - CASE 2

syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda t
 
p_H0_T = p_H0 * (1 + t);
p_L0_T = p_L0 * (1 + t);
 
U_H1 = R + theta*s_H - p_H1 - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0_T - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0_T - tau*x^2;
 
disp(U_H1);
disp(U_H0);
disp(U_L1);
disp(U_L0);

Probabilities - CASE 2

denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
 
disp(Pr_H1)
disp(Pr_H0)
disp(Pr_L1)
disp(Pr_L0)
 
% Pr_H1 = exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))));
% Pr_H0 = exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))));
% Pr_L1 = exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))));
% Pr_L0 = exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))));

Demand functions - CASE 2

syms x theta p_H1 p_H0 p_L1 p_L0 theta_min theta_max
 
f_x = 1;
g_theta = 1 / (theta_max - theta_min);
 
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
 
disp('D_H1 = '); disp(D_H1)
D_H1 =
disp('D_H0 = '); disp(D_H0)
D_H0 =
disp('D_L1 = '); disp(D_L1)
D_L1 =
disp('D_L0 = '); disp(D_L0)
D_L0 =
 
% D_H1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% D_H0 = int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% D_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% D_L0 = int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);

Profit functions - CASE 2

Pi_ij = (p_ij - c_ij) * D_ij
 
% tax just applied in demand!
 
syms c_H1 c_H0 c_L1 c_L0
 
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
 
disp('Pi_H1 = '); disp(Pi_H1)
Pi_H1 =
disp('Pi_H0 = '); disp(Pi_H0)
Pi_H0 =
disp('Pi_L1 = '); disp(Pi_L1)
Pi_L1 =
disp('Pi_L0 = '); disp(Pi_L0)
Pi_L0 =
 
% Pi_H1 = -int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_H1 - p_H1);
% Pi_H0 = -(c_H0 - p_H0)*int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% Pi_L1 = -int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_L1 - p_L1);
% Pi_L0 = -(c_L0 - p_L0)*int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);

Profit maximization - CASE 2

d(Pi_ij) / d(p_ij) = 0
dPi_H1_dpH1 = diff(Pi_H1, p_H1) == 0;
dPi_H0_dpH0 = diff(Pi_H0, p_H0) == 0;
dPi_L1_dpL1 = diff(Pi_L1, p_L1) == 0;
dPi_L0_dpL0 = diff(Pi_L0, p_L0) == 0;
 
disp('dPi_H1_dpH1:'); disp(dPi_H1_dpH1)
dPi_H1_dpH1:
disp('dPi_H0_dpH0:'); disp(dPi_H0_dpH0)
dPi_H0_dpH0:
disp('dPi_L1_dpL1:'); disp(dPi_L1_dpL1)
dPi_L1_dpL1:
disp('dPi_L0_dpL0:'); disp(dPi_L0_dpL0)
dPi_L0_dpL0:
 
% dPi_H1_dpH1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_H1 - p_H1)*int(int((lambda*exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))) - (lambda*exp(2*lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
% dPi_H0_dpH0 = int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - int(int((lambda*exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))) - (lambda*exp(2*lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1)*(c_H0 - p_H0) == 0;
% dPi_L1_dpL1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_L1 - p_L1)*int(int((lambda*exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))) - (lambda*exp(2*lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
% dPi_L0_dpL0 = int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - int(int((lambda*exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))) - (lambda*exp(2*lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1)*(c_L0 - p_L0) == 0;

Solve FOC - CASE 2

syms p_H1 p_H0 p_L1 p_L0 theta x
 
R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
theta_min = 0.2;
theta_max = 2.2;
t = 0.25;
c_H1 = 0.35;
c_H0 = 0.35;
c_L1 = 0.1;
c_L0 = 0.1;
 
denominator = (theta_max - theta_min) * ( ...
exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + ...
exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + ...
exp(lambda*(-tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + ...
exp(lambda*(-tau*x^2 + R + s_L*theta - p_L0*(t + 1))) );
 
% FOCs
foc_H1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_H1 - p_H1) * int(int((lambda * exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_H1 - tau*(x - 1)^2 + s_H*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
 
foc_H0 = int(int(-exp(lambda*(-tau*x^2 + R + s_H*theta - p_H0*(t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- int(int((lambda * exp(lambda*(-tau*x^2 + R + s_H*theta - p_H0*(t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (-tau*x^2 + R + s_H*theta - p_H0*(t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) * (c_H0 - p_H0) == 0;
 
foc_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_L1 - p_L1) * int(int((lambda * exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_L1 - tau*(x - 1)^2 + s_L*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
 
foc_L0 = int(int(-exp(lambda*(-tau*x^2 + R + s_L*theta - p_L0*(t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- int(int((lambda * exp(lambda*(-tau*x^2 + R + s_L*theta - p_L0*(t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (-tau*x^2 + R + s_L*theta - p_L0*(t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) * (c_L0 - p_L0) == 0;
 
solution = vpasolve([foc_H1, foc_H0, foc_L1, foc_L0], [p_H1, p_H0, p_L1, p_L0]);
 
disp('Optimal Prices:');
Optimal Prices:
disp(['p_H1** = ', char(solution.p_H1)]);
p_H1** = 0.98245738528891675184015215351829
disp(['p_H0** = ', char(solution.p_H0)]);
p_H0** = 0.8451818220778982201818734232353
disp(['p_L1** = ', char(solution.p_L1)]);
p_L1** = 0.67507753269566794239079645081273
disp(['p_L0** = ', char(solution.p_L0)]);
p_L0** = 0.5602969265225777050338138212517

Numerical results - CASE 2

R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
theta_min = 0.2;
theta_max = 2.2;
c_H1 = 0.35;
c_H0 = 0.35;
c_L1 = 0.1;
c_L0 = 0.1;
t = 0.25;
 
p_H1 = 0.9824;
p_H0 = 0.8451;
p_L1 = 0.6750;
p_L0 = 0.5602;
 
denominator = @(x, theta) (theta_max - theta_min) .* ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0*(1+t))) + ...
exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0*(1+t))) );
 
D_H1 = integral2(@(theta, x) exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0*(1+t))) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0*(1+t))) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
 
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
 
fprintf('Demands:\n');
Demands:
fprintf('D_H1 = %.4f\n', D_H1);
D_H1 = 0.3298
fprintf('D_H0 = %.4f\n', D_H0);
D_H0 = 0.2923
fprintf('D_L1 = %.4f\n', D_L1);
D_L1 = 0.1913
fprintf('D_L0 = %.4f\n', D_L0);
D_L0 = 0.1866
 
fprintf('\nProfits:\n');
Profits:
fprintf('Pi_H1 = %.4f\n', Pi_H1);
Pi_H1 = 0.2086
fprintf('Pi_H0 = %.4f\n', Pi_H0);
Pi_H0 = 0.1447
fprintf('Pi_L1 = %.4f\n', Pi_L1);
Pi_L1 = 0.1100
fprintf('Pi_L0 = %.4f\n', Pi_L0);
Pi_L0 = 0.0859
 
p_H0_tax = p_H0 * (1 + t);
p_L0_tax = p_L0 * (1 + t);
 
fprintf('\nConsumer Prices:\n');
Consumer Prices:
fprintf('p_H0_VAT = %.4f\n', p_H0_tax);
p_H0_VAT = 1.0564
fprintf('p_L0_VAT = %.4f\n', p_L0_tax);
p_L0_VAT = 0.7003
 
 
% end of CASE 2:
clc; clear; close all;

Master-Thesis - Model adjusted for VAT on ICEV and high qualiy EV:

Utility - CASE 3

syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda t
 
p_H0_T = p_H0 * (1 + t);
p_L0_T = p_L0 * (1 + t);
p_H1_T = p_H1 * (1 + t);
 
 
U_H1 = R + theta*s_H - p_H1_T - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0_T - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0_T - tau*x^2;
 
disp(U_H1);
disp(U_H0);
disp(U_L1);
disp(U_L0);

Probabilities - CASE 3

denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
 
disp(Pr_H1)
disp(Pr_H0)
disp(Pr_L1)
disp(Pr_L0)

Demand functions - CASE 3

syms x theta p_H1 p_H0 p_L1 p_L0 theta_min theta_max
 
f_x = 1;
g_theta = 1 / (theta_max - theta_min);
 
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
 
disp('D_H1 = '); disp(D_H1)
D_H1 =
disp('D_H0 = '); disp(D_H0)
D_H0 =
disp('D_L1 = '); disp(D_L1)
D_L1 =
disp('D_L0 = '); disp(D_L0)
D_L0 =
 
% D_H1 = int(int(-exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1
% D_H0 = int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)
% D_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)
% D_L0 = int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)

Profit functions - CASE 3

Pi_ij = (p_ij - c_ij) * D_ij
 
syms c_H1 c_H0 c_L1 c_L0
 
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
 
disp('Pi_H1 = '); disp(Pi_H1)
Pi_H1 =
disp('Pi_H0 = '); disp(Pi_H0)
Pi_H0 =
disp('Pi_L1 = '); disp(Pi_L1)
Pi_L1 =
disp('Pi_L0 = '); disp(Pi_L0)
Pi_L0 =
 
% Pi_H1 = -(c_H1 - p_H1)*int(int(-exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% Pi_H0 = -int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_H0 - p_H0);
% Pi_L1 = -int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_L1 - p_L1);
% Pi_L0 = -int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_L0 - p_L0);

Profit maximization - CASE 3

d(Pi_ij) / d(p_ij) = 0
 
dPi_H1_dpH1 = diff(Pi_H1, p_H1) == 0;
dPi_H0_dpH0 = diff(Pi_H0, p_H0) == 0;
dPi_L1_dpL1 = diff(Pi_L1, p_L1) == 0;
dPi_L0_dpL0 = diff(Pi_L0, p_L0) == 0;
 
disp('dPi_H1_dpH1 = '); disp(dPi_H1_dpH1)
dPi_H1_dpH1 =
disp('dPi_H0_dpH0 = '); disp(dPi_H0_dpH0)
dPi_H0_dpH0 =
disp('dPi_L1_dpL1 = '); disp(dPi_L1_dpL1)
dPi_L1_dpL1 =
disp('dPi_L0_dpL0 = '); disp(dPi_L0_dpL0)
dPi_L0_dpL0 =
 
% dPi_H1_dpH1 = int(int(-exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_H1 - p_H1)*int(int((lambda*exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))) - (lambda*exp(2*lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
% dPi_H0_dpH0 = int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_H0 - p_H0)*int(int((lambda*exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))) - (lambda*exp(2*lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
% dPi_L1_dpL1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - int(int((lambda*exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))) - (lambda*exp(2*lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1)*(c_L1 - p_L1) == 0;
% dPi_L0_dpL0 = int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_L0 - p_L0)*int(int((lambda*exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))) - (lambda*exp(2*lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;

Solve FOC - CASE 3

syms p_H1 p_H0 p_L1 p_L0 theta x
 
R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
theta_min = 0.2;
theta_max = 2.2;
t = 0.25;
c_H1 = 0.35;
c_H0 = 0.35;
c_L1 = 0.1;
c_L0 = 0.1;
 
denominator = (theta_max - theta_min) * ( ...
exp(lambda * (R - p_L1 - tau * (x - 1)^2 + s_L * theta)) + ...
exp(lambda * (-tau * x^2 + R + s_H * theta - p_H0 * (t + 1))) + ...
exp(lambda * (-tau * x^2 + R + s_L * theta - p_L0 * (t + 1))) + ...
exp(lambda * (R - tau * (x - 1)^2 + s_H * theta - p_H1 * (t + 1))) );
 
% FOCs:
foc_H1 = int(int(-exp(lambda * (R - tau * (x - 1)^2 + s_H * theta - p_H1 * (t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_H1 - p_H1) * int(int((lambda * exp(lambda * (R - tau * (x - 1)^2 + s_H * theta - p_H1 * (t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (R - tau * (x - 1)^2 + s_H * theta - p_H1 * (t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
 
foc_H0 = int(int(-exp(lambda * (-tau * x^2 + R + s_H * theta - p_H0 * (t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_H0 - p_H0) * int(int((lambda * exp(lambda * (-tau * x^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (-tau * x^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
 
foc_L1 = int(int(-exp(lambda * (R - p_L1 - tau * (x - 1)^2 + s_L * theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_L1 - p_L1) * int(int((lambda * exp(lambda * (R - p_L1 - tau * (x - 1)^2 + s_L * theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_L1 - tau * (x - 1)^2 + s_L * theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
 
foc_L0 = int(int(-exp(lambda * (-tau * x^2 + R + s_L * theta - p_L0 * (t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_L0 - p_L0) * int(int((lambda * exp(lambda * (-tau * x^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (-tau * x^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
 
solution = vpasolve([foc_H1, foc_H0, foc_L1, foc_L0], [p_H1, p_H0, p_L1, p_L0]);
 
disp('Optimal Prices:');
Optimal Prices:
disp(['p_H1*** = ', char(solution.p_H1)]);
p_H1*** = 0.84730365281426619546342786138934
disp(['p_H0*** = ', char(solution.p_H0)]);
p_H0*** = 0.84809017803927929379959703727505
disp(['p_L1*** = ', char(solution.p_L1)]);
p_L1*** = 0.67999498272950621551040317326515
disp(['p_L0*** = ', char(solution.p_L0)]);
p_L0*** = 0.56166362780346060630817744594582
 

Numerical results - CASE 3

R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
c_H1 = 0.35;
c_H0 = 0.35;
c_L1 = 0.1;
c_L0 = 0.1;
theta_min = 0.2;
theta_max = 2.2;
t = 0.25;
 
p_H1 = 0.8473;
p_H0 = 0.8480;
p_L1 = 0.6799;
p_L0 = 0.5616;
 
p_H1_tax = p_H1 * (1 + t);
p_H0_tax = p_H0 * (1 + t);
p_L0_tax = p_L0 * (1 + t);
 
denominator = @(x, theta) (theta_max - theta_min) .* ( ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0_tax)) + ...
exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0_tax)) + ...
exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1_tax)));
 
D_H1 = integral2(@(theta, x) exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1_tax)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0_tax)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0_tax)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
 
fprintf('Demands:\n');
Demands:
fprintf('D_H1 = %.4f\n', D_H1);
D_H1 = 0.3023
fprintf('D_H0 = %.4f\n', D_H0);
D_H0 = 0.3031
fprintf('D_L1 = %.4f\n', D_L1);
D_L1 = 0.2012
fprintf('D_L0 = %.4f\n', D_L0);
D_L0 = 0.1934
 
fprintf('\nProfits:\n');
Profits:
fprintf('Pi_H1 = %.4f\n', Pi_H1);
Pi_H1 = 0.1504
fprintf('Pi_H0 = %.4f\n', Pi_H0);
Pi_H0 = 0.1509
fprintf('Pi_L1 = %.4f\n', Pi_L1);
Pi_L1 = 0.1167
fprintf('Pi_L0 = %.4f\n', Pi_L0);
Pi_L0 = 0.0893
 
fprintf('\nConsumer Prices:\n');
Consumer Prices:
fprintf('p_H1_VAT = %.4f\n', p_H1_tax);
p_H1_VAT = 1.0591
fprintf('p_H0_VAT = %.4f\n', p_H0_tax);
p_H0_VAT = 1.0600
fprintf('p_L0_VAT = %.4f\n', p_L0_tax);
p_L0_VAT = 0.7020
 

Master-Thesis - Model calculations for top and bottom 10% WTP

clc; clear; close all;

Numerical results

CASE 1
R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
p_H1 = 0.9781;
p_H0 = 0.9781;
p_L1 = 0.6729;
p_L0 = 0.6729;
 
theta_min = 0.2;
theta_max = 2.2;
theta_90 = theta_min + 0.9 * (theta_max - theta_min);
theta_10 = theta_min + 0.1 * (theta_max - theta_min);
 
denominator = @(theta,x) (exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) + exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) );
 
D_H1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta,x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) ./ denominator(theta,x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta,x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) ./ denominator(theta,x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta,x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0)) ./ denominator(theta,x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta,x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0)) ./ denominator(theta,x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
fprintf('D_H1 (Top 10%% WTP) = %.4f\n', D_H1_top10);
D_H1 (Top 10% WTP) = 0.2040
fprintf('D_H0 (Top 10%% WTP) = %.4f\n', D_H0_top10);
D_H0 (Top 10% WTP) = 0.2040
fprintf('D_L1 (Top 10%% WTP) = %.4f\n', D_L1_top10);
D_L1 (Top 10% WTP) = 0.0460
fprintf('D_L0 (Top 10%% WTP) = %.4f\n', D_L0_top10);
D_L0 (Top 10% WTP) = 0.0460
 
fprintf('D_H1 (Bottom 10%% WTP) = %.4f\n', D_H1_low10);
D_H1 (Bottom 10% WTP) = 0.1058
fprintf('D_H0 (Bottom 10%% WTP) = %.4f\n', D_H0_low10);
D_H0 (Bottom 10% WTP) = 0.1058
fprintf('D_L1 (Bottom 10%% WTP) = %.4f\n', D_L1_low10);
D_L1 (Bottom 10% WTP) = 0.1442
fprintf('D_L0 (Bottom 10%% WTP) = %.4f\n', D_L0_low10);
D_L0 (Bottom 10% WTP) = 0.1442

CASE 2

R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
t = 0.25;
p_H1 = 0.9824;
p_H0 = 0.8451;
p_L1 = 0.6750;
p_L0 = 0.5602;
 
theta_min = 0.2;
theta_max = 2.2;
theta_90 = theta_min + 0.9 * (theta_max - theta_min);
theta_10 = theta_min + 0.1 * (theta_max - theta_min);
 
 
denominator = @(theta, x) ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) + ...
exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) );
 
D_H1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
fprintf('D_H1 (Top 10%% WTP) = %.4f\n', D_H1_top10);
D_H1 (Top 10% WTP) = 0.2138
fprintf('D_H0 (Top 10%% WTP) = %.4f\n', D_H0_top10);
D_H0 (Top 10% WTP) = 0.1902
fprintf('D_L1 (Top 10%% WTP) = %.4f\n', D_L1_top10);
D_L1 (Top 10% WTP) = 0.0485
fprintf('D_L0 (Top 10%% WTP) = %.4f\n', D_L0_top10);
D_L0 (Top 10% WTP) = 0.0475
 
fprintf('D_H1 (Bottom 10%% WTP) = %.4f\n', D_H1_low10);
D_H1 (Bottom 10% WTP) = 0.1091
fprintf('D_H0 (Bottom 10%% WTP) = %.4f\n', D_H0_low10);
D_H0 (Bottom 10% WTP) = 0.0962
fprintf('D_L1 (Bottom 10%% WTP) = %.4f\n', D_L1_low10);
D_L1 (Bottom 10% WTP) = 0.1495
fprintf('D_L0 (Bottom 10%% WTP) = %.4f\n', D_L0_low10);
D_L0 (Bottom 10% WTP) = 0.1452

CASE 3

R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
t = 0.25;
p_H1 = 0.8473;
p_H0 = 0.8480;
p_L1 = 0.6799;
p_L0 = 0.5616;
 
theta_min = 0.2;
theta_max = 2.2;
theta_90 = theta_min + 0.9 * (theta_max - theta_min);
theta_10 = theta_min + 0.1 * (theta_max - theta_min);
 
denominator = @(theta, x) (exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1 * (1 + t))) + exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) + exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) );
 
D_H1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 * (1 + t) - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 * (1 + t) - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
 
fprintf('D_H1 (Top 10%% WTP) = %.4f\n', D_H1_top10);
D_H1 (Top 10% WTP) = 0.1989
fprintf('D_H0 (Top 10%% WTP) = %.4f\n', D_H0_top10);
D_H0 (Top 10% WTP) = 0.1991
fprintf('D_L1 (Top 10%% WTP) = %.4f\n', D_L1_top10);
D_L1 (Top 10% WTP) = 0.0521
fprintf('D_L0 (Top 10%% WTP) = %.4f\n', D_L0_top10);
D_L0 (Top 10% WTP) = 0.0499
 
fprintf('D_H1 (Bottom 10%% WTP) = %.4f\n', D_H1_low10);
D_H1 (Bottom 10% WTP) = 0.0979
fprintf('D_H0 (Bottom 10%% WTP) = %.4f\n', D_H0_low10);
D_H0 (Bottom 10% WTP) = 0.0983
fprintf('D_L1 (Bottom 10%% WTP) = %.4f\n', D_L1_low10);
D_L1 (Bottom 10% WTP) = 0.1547
fprintf('D_L0 (Bottom 10%% WTP) = %.4f\n', D_L0_low10);
D_L0 (Bottom 10% WTP) = 0.1490
 

clc; clear; close all;

Master-Thesis - Model calculations for the elasticities

Elasticities CASE 1

syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda theta_min theta_max
 
U_H1 = R + theta*s_H - p_H1 - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0 - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0 - tau*x^2;
 
disp(U_H1);
disp(U_H0);
disp(U_L1);
disp(U_L0);
 
denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
 
disp(Pr_H1)
disp(Pr_H0)
disp(Pr_L1)
disp(Pr_L0)
 
f_x = 1;
g_theta = 1 / (theta_max - theta_min);
 
 
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
 
disp('D_H1 = '); disp(D_H1)
D_H1 =
disp('D_H0 = '); disp(D_H0)
D_H0 =
disp('D_L1 = '); disp(D_L1)
D_L1 =
disp('D_L0 = '); disp(D_L0)
D_L0 =
 
dD_H1_dp_H1 = diff(D_H1, p_H1);
dD_H0_dp_H0 = diff(D_H0, p_H0);
dD_L1_dp_L1 = diff(D_L1, p_L1);
dD_L0_dp_L0 = diff(D_L0, p_L0);
 
disp(dD_H1_dp_H1)
disp(dD_H0_dp_H0)
disp(dD_L1_dp_L1)
disp(dD_L0_dp_L0)
 
e_H1 = simplify(dD_H1_dp_H1 * (p_H1/D_H1));
e_H0 = simplify(dD_H0_dp_H0 * (p_H0/D_H0));
e_L1 = simplify(dD_L1_dp_L1 * (p_L1/D_L1));
e_L0 = simplify(dD_L0_dp_L0 * (p_L0/D_L0));
 
disp('e_H1 ='); disp(e_H1)
e_H1 =
disp('e_H0 ='); disp(e_H0)
e_H0 =
disp('e_L1 ='); disp(e_L1)
e_L1 =
disp('e_L0 ='); disp(e_L0)
e_L0 =

Nummerical calculation elasticity CASE 1

R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
theta_min = 0.2;
theta_max = 2.2;
p_H1 = 0.9781;
p_H0 = 0.9781;
p_L1 = 0.6729;
p_L0 = 0.6729;
 
 
denominator = @(x, theta) ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) + ...
exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) );
 
D_H1 = integral2(@(theta, x) - exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_H1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_H0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R - p_H0 + s_H * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_L1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_L0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R - p_L0 + s_L * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
 
e_H1 = dD_H1_dP * (p_H1 / D_H1);
e_H0 = dD_H0_dP * (p_H0 / D_H0);
e_L1 = dD_L1_dP * (p_L1 / D_L1);
e_L0 = dD_L0_dP * (p_L0 / D_L0);
 
 
disp('Computed Demand Values:')
Computed Demand Values:
fprintf('D_H1 = %.6f\n', D_H1);
D_H1 = 0.316600
fprintf('D_H0 = %.6f\n', D_H0);
D_H0 = 0.316600
fprintf('D_L1 = %.6f\n', D_L1);
D_L1 = 0.183400
fprintf('D_L0 = %.6f\n', D_L0);
D_L0 = 0.183400
 
disp('Computed Demand Derivatives:')
Computed Demand Derivatives:
fprintf('dD_H1/dP = %.6f\n', dD_H1_dP);
dD_H1/dP = -0.762406
fprintf('dD_H0/dP = %.6f\n', dD_H0_dP);
dD_H0/dP = -0.762406
fprintf('dD_L1/dP = %.6f\n', dD_L1_dP);
dD_L1/dP = -0.413527
fprintf('dD_L0/dP = %.6f\n', dD_L0_dP);
dD_L0/dP = -0.413527
 
disp('Computed Elasticities:')
Computed Elasticities:
fprintf('e_H1 = %.4f\n', e_H1);
e_H1 = -2.3554
fprintf('e_H0 = %.4f\n', e_H0);
e_H0 = -2.3554
fprintf('e_L1 = %.4f\n', e_L1);
e_L1 = -1.5172
fprintf('e_L0 = %.4f\n', e_L0);
e_L0 = -1.5172
 

CASE 2

clc; clear; close all;
 
syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda theta_min theta_max t
 
U_H1 = R + theta*s_H - p_H1 - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0*(1+t) - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0*(1+t) - tau*x^2;
 
disp(U_H1);
disp(U_H0);
disp(U_L1);
disp(U_L0);
 
denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
 
disp(Pr_H1)
disp(Pr_H0)
disp(Pr_L1)
disp(Pr_L0)
 
f_x = 1;
g_theta = 1 / (theta_max - theta_min);
 
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
 
disp('D_H1:')
D_H1:
disp(D_H1)
disp('D_H0:')
D_H0:
disp(D_H0)
disp('D_L1:')
D_L1:
disp(D_L1)
disp('D_L0:')
D_L0:
disp(D_L0)
 
dD_H1_dp_H1 = diff(D_H1, p_H1);
dD_H0_dp_H0 = diff(D_H0, p_H0);
dD_L1_dp_L1 = diff(D_L1, p_L1);
dD_L0_dp_L0 = diff(D_L0, p_L0);
 
disp(dD_H1_dp_H1)
disp(dD_H0_dp_H0)
disp(dD_L1_dp_L1)
disp(dD_L0_dp_L0)
 
e_H1 = simplify(dD_H1_dp_H1 * (p_H1/D_H1));
e_H0 = simplify(dD_H0_dp_H0 * (p_H0/D_H0));
e_L1 = simplify(dD_L1_dp_L1 * (p_L1/D_L1));
e_L0 = simplify(dD_L0_dp_L0 * (p_L0/D_L0));
 
disp('e_H1 ='); disp(e_H1)
e_H1 =
disp('e_H0 ='); disp(e_H0)
e_H0 =
disp('e_L1 ='); disp(e_L1)
e_L1 =
disp('e_L0 ='); disp(e_L0)
e_L0 =

Nummerical calculation elasticity CASE 2

 
R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
theta_min = 0.2;
theta_max = 2.2;
t = 0.25;
p_H1 = 0.9824;
p_H0 = 0.8451;
p_L1 = 0.6750;
p_L0 = 0.5616;
 
 
denominator = @(x, theta) ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (R - tau * x.^2 + s_H * theta - (1 + t) * p_H0)) + ...
exp(lambda * (R - tau * x.^2 + s_L * theta - (1 + t) * p_L0)) );
 
 
D_H1 = integral2(@(theta, x) exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0 = integral2(@(theta, x) exp(lambda * (R - tau * x.^2 + s_H * theta - (1 + t) * p_H0)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0 = integral2(@(theta, x) exp(lambda * (R - tau * x.^2 + s_L * theta - (1 + t) * p_L0)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_H1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_H0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_L1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_L0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
e_H1 = dD_H1_dP * (p_H1 / D_H1);
e_H0 = dD_H0_dP * (p_H0 / D_H0);
e_L1 = dD_L1_dP * (p_L1 / D_L1);
e_L0 = dD_L0_dP * (p_L0 / D_L0);
 
disp('Computed Demand Values with Tax:')
Computed Demand Values with Tax:
fprintf('D_H1 = %.6f\n', D_H1);
D_H1 = 0.329951
fprintf('D_H0 = %.6f\n', D_H0);
D_H0 = 0.292477
fprintf('D_L1 = %.6f\n', D_L1);
D_L1 = 0.191448
fprintf('D_L0 = %.6f\n', D_L0);
D_L0 = 0.186124
 
disp('Computed Demand Derivatives with Tax:')
Computed Demand Derivatives with Tax:
fprintf('dD_H1/dP = %.6f\n', dD_H1_dP);
dD_H1/dP = -0.798134
fprintf('dD_H0/dP = %.6f\n', dD_H0_dP);
dD_H0/dP = -0.871860
fprintf('dD_L1/dP = %.6f\n', dD_L1_dP);
dD_L1/dP = -0.432904
fprintf('dD_L0/dP = %.6f\n', dD_L0_dP);
dD_L0/dP = -0.526133
 
disp('Computed Elasticities with Tax:')
Computed Elasticities with Tax:
fprintf('e_H1 = %.4f\n', e_H1);
e_H1 = -2.3764
fprintf('e_H0 = %.4f\n', e_H0);
e_H0 = -2.5192
fprintf('e_L1 = %.4f\n', e_L1);
e_L1 = -1.5263
fprintf('e_L0 = %.4f\n', e_L0);
e_L0 = -1.5875
 

CASE 3

clc; clear; close all;
 
syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda theta_min theta_max t
 
U_H1 = R + theta*s_H - p_H1*(1+t) - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0*(1+t) - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0*(1+t) - tau*x^2;
 
disp(U_H1);
disp(U_H0);
disp(U_L1);
disp(U_L0);
 
 
denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
 
disp(Pr_H1)
disp(Pr_H0)
disp(Pr_L1)
disp(Pr_L0)
 
 
f_x = 1;
g_theta = 1 / (theta_max - theta_min);
 
 
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
 
disp('D_H1 = '); disp(D_H1)
D_H1 =
disp('D_H0 = '); disp(D_H0)
D_H0 =
disp('D_L1 = '); disp(D_L1)
D_L1 =
disp('D_L0 = '); disp(D_L0)
D_L0 =
 
dD_H1_dp_H1 = diff(D_H1, p_H1);
dD_H0_dp_H0 = diff(D_H0, p_H0);
dD_L1_dp_L1 = diff(D_L1, p_L1);
dD_L0_dp_L0 = diff(D_L0, p_L0);
 
disp(dD_H1_dp_H1)
disp(dD_H0_dp_H0)
disp(dD_L1_dp_L1)
disp(dD_L0_dp_L0)
 
e_H1 = simplify(dD_H1_dp_H1 * (p_H1/D_H1));
e_H0 = simplify(dD_H0_dp_H0 * (p_H0/D_H0));
e_L1 = simplify(dD_L1_dp_L1 * (p_L1/D_L1));
e_L0 = simplify(dD_L0_dp_L0 * (p_L0/D_L0));
 
 
disp('e_H1 ='); disp(e_H1)
e_H1 =
disp('e_H0 ='); disp(e_H0)
e_H0 =
disp('e_L1 ='); disp(e_L1)
e_L1 =
disp('e_L0 ='); disp(e_L0)
e_L0 =

Numerical calculations elasticity CASE 3

R = 1;
lambda = 2;
tau = 1;
s_H = 0.7;
s_L = 0.2;
theta_min = 0.2;
theta_max = 2.2;
t = 0.25;
p_H1 = 0.8473;
p_H0 = 0.8480;
p_L1 = 0.6799;
p_L0 = 0.5616;
 
denominator = @(x, theta) ( ...
exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - (1 + t) * p_H1)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (R - tau * x.^2 + s_H * theta - (1 + t) * p_H0)) + ...
exp(lambda * (R - tau * x.^2 + s_L * theta - (1 + t) * p_L0)) );
 
D_H1 = integral2(@(theta, x) ...
-exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1 * (t + 1))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_H0 = integral2(@(theta, x) ...
-exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L1 = integral2(@(theta, x) ...
-exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
D_L0 = integral2(@(theta, x) ...
-exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
 
dD_H1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_H0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_L1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
dD_L0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
 
 
e_H1 = dD_H1_dP * (p_H1 / D_H1);
e_H0 = dD_H0_dP * (p_H0 / D_H0);
e_L1 = dD_L1_dP * (p_L1 / D_L1);
e_L0 = dD_L0_dP * (p_L0 / D_L0);
 
 
disp('Computed Demand Values with Tax:')
Computed Demand Values with Tax:
fprintf('D_H1 = %.6f\n', D_H1);
D_H1 = 0.302344
fprintf('D_H0 = %.6f\n', D_H0);
D_H0 = 0.303062
fprintf('D_L1 = %.6f\n', D_L1);
D_L1 = 0.201178
fprintf('D_L0 = %.6f\n', D_L0);
D_L0 = 0.193416
 
disp('Computed Demand Derivatives with Tax:')
Computed Demand Derivatives with Tax:
fprintf('dD_H1/dP = %.6f\n', dD_H1_dP);
dD_H1/dP = -0.903740
fprintf('dD_H0/dP = %.6f\n', dD_H0_dP);
dD_H0/dP = -0.906869
fprintf('dD_L1/dP = %.6f\n', dD_L1_dP);
dD_L1/dP = -0.457854
fprintf('dD_L0/dP = %.6f\n', dD_L0_dP);
dD_L0/dP = -0.548126
 
disp('Computed Elasticities with Tax:')
Computed Elasticities with Tax:
fprintf('e_H1 = %.4f\n', e_H1);
e_H1 = -2.5327
fprintf('e_H0 = %.4f\n', e_H0);
e_H0 = -2.5375
fprintf('e_L1 = %.4f\n', e_L1);
e_L1 = -1.5474
fprintf('e_L0 = %.4f\n', e_L0);
e_L0 = -1.5915